Unveiling the Emergent Traits of Chiral Spin Textures in Magnetic Multilayers

Abstract Magnetic skyrmions are topologically wound nanoscale textures of spins whose ambient stability and electrical manipulation in multilayer films have led to an explosion of research activities. While past efforts focused predominantly on isolated skyrmions, recently ensembles of chiral spin textures, consisting of skyrmions and magnetic stripes, are shown to possess rich interactions with potential for device applications. However, several fundamental aspects of chiral spin texture phenomenology remain to be elucidated, including their domain wall (DW) structure, thermodynamic stability, and morphological transitions. Here the evolution of these textural characteristics are unveiled on a tunable multilayer platform—wherein chiral interactions governing spin texture energetics can be widely varied—using a combination of full‐field electron and soft X‐ray microscopies with numerical simulations. With increasing chiral interactions, the emergence of Néel helicity, followed by a marked reduction in domain compressibility, and finally a transformation in the skyrmion formation mechanism are demonstrated. Together with an analytical model, these experiments establish a comprehensive microscopic framework for investigating and tailoring chiral spin texture character in multilayer films.


S1. Magnetic Properties of Multilayer Samples
Magnetic Properties. The magnetic properties of the samples studied in this work are tabulated in Table S1 and plotted in Figure S1. The saturation magnetization, M s , and effective anisotropy, K eff , were determined from magnetometry measurements (see Figure S2a-d). Following procedures established in previous works 1 , the exchange stiffness, A est , was estimated with density functional theory (DFT) 2 . Meanwhile, the iDMI value, D est was estimated by comparing the periodicity of micromagnetic simulations at zero field with that from magnetic force microscopy (MFM) images (see Figure S2q-t) 1,3,4 . Notably, the values A est and D est for the 14x multilayers studied in this work are in good agreement with values obtained for corresponding 1x films from Brillouin light scattering (BLS) experiments 2 . To check if a small iDMI affects the simulation of the S Co(10) stack, we performed additional hysteresis loop simulations with iDMI of 0.1 mJ/m 22 and found the results to be virtually identical to our original zero iDMI simulations.
Magnetometry and magnetic microscopy. In Figure S2, we show the magnetization hysteresis loops, M (H), with applied field µ 0 H in the in-plane (IP) and out-of-plane (OP) directions. We also show an overview of LTEM, MTXM and micromagnetic simulations for the four samples at zero field or low fields, in the labyrinthine stripe state. The  Table S1. Magnetic properties of multilayer samples. List of multilayer samples used in this work, with layer thickness in angstroms in parentheses. The magnetic properties -saturation magnetization (Ms), effective magnetic anisotropy (K eff ), estimated iDMI (Dest), estimated exchange stiffness (Aest) as well as the stability parameter, κ -are tabulated. magnetic morphology captured by the two experimental imaging techniques is consistent. The simulation results are also largely in agreement with the experimental images.

S2. Image analysis
Basic image processing. The analysis of LTEM and MTXM images used in this work was performed with customwritten Python code. To remove low spatial frequency background, a duplicate convolved with a large Gaussian kernel was subtracted from the original image. To reduce the high spatial frequency noise while preserving as much information as possible, a small median filter was then applied. In cases where binarization was required, the threshold was automatically selected using the Otsu algorithm. To remove further noise, if any, small regions below a certain areal threshold were disregarded.
Domain width analysis. Manuscript Fig. 2 shows the averaged domain width from MTXM images. To obtain this result, background removal, denoising and binarization were performed as described above. Next, a Euclidean distance transform -which yields the shortest distance from the domain edge to the medial axis -was applied on the binarized image using built-in functions in the scikit-image library 5 . Finally, the average domain width was extracted as twice the average distance to the medial axis.
Linecut analysis. Manuscript Fig. 1 shows linecuts across stripe domains at fixed intervals from LTEM images. To obtain these linecuts, a skeletonization operation was enacted on the LTEM images after the basic image processing and binarization described above. Next, all foreground pixels connected to more than two other pixels were removed. This ensured that there were no branches and that all regions, defined as foreground pixels that are 4-connected, are linear. Each region was then fitted to a spline, and the normal to the spline was computed.
Finally, linecuts were acquired from the original, unprocessed LTEM image using methods in scikit-image library, at regular intervals on the spline. Linecuts exceeding the boundaries of the image were excluded. The position and the direction of the algorithmically chosen linecuts are shown in Figure S3. The acquired linecuts were further binned by angle, and only those within 10°from the tilt axis were included in calculating the averaged output. For analysis of LTEM images acquired at varying tilt angles, an additional preprocessing step was required. The sequence of images had to be aligned to the first image, as the latter was used to determine the position of each linecut. The image alignment was done manually with ImageJ by selecting about 15 corresponding landmarks on each image and calculating an affine transformation to map each image in the tilt sequence to the first image.
Manuscript Fig. 1 does not show the linecut for zero tilt of Fe(0)/Co (10). This is because, in the absence of visible magnetic textures at zero tilt, the alignment landmarks for that image could not be located. Similar image sequence alignment issues precluded the use of this method for analyzing Fe(3)/Co (7) images (see § S3).   identified as above is shown in Figure S4.

S3. LTEM Imaging
LTEM contrast simulation. In this work, spatial evolution of LTEM contrast of textures has been used to quantify their helicity. As LTEM contrast is not intuitive to non-experts, we show in Figure S5 the simulated LTEM contrast -generated using an open-source software, MALTS 6 -for an idealized magnetic stripe with Bloch and Néel DWs.
Also shown are the linecuts perpendicular to the domain for various tilt angles (c.f. manuscript Fig. 1f-g). At zero tilt, the Bloch DW pair produces contrast symmetric about the center of the stripe, whereas the Néel DW pair give zero contrast. In both cases, an additional antisymmetric contrast is present for finite tilt.
Additional LTEM data for DW helicity. In manuscript Fig. 1, we have shown tilt-dependent LTEM data for S Co(10), confirming Bloch helicity of its DWs, and for Fe(0)/Co(10) -wherein DWs have Néel helicity. In Figure S6, we show for comparison LTEM data for Fe(2)/Co (8), which also has, unsurprisingly, Néel DWs. Finally, the linecut extraction technique used for the previous three samples could not be applied to Fe(3)/Co (7). In this case, owing to the higher texture density, images with varying tilt could not be reliably aligned (see § S2). Nonetheless, a visual inspection of the contrast evolution with tilt suggests that the textural helicity for Fe(3)/Co (7) is very similar to that of Fe(2)/Co (8) hosting Néel DWs.  Transport-of-Intensity Equation (TIE). The transport-of-intensity equation (TIE) has previously been used to reconstruct the magnetization of Bloch and Néel textures with varying degrees of success 7 . Correspondingly in Figure S7, we show the IP magnetic induction for S Co(10) and Fe(0)/Co(10), each reconstructed using two LTEM images at defoci of ±2 mm using a TIE-based commercial software, QPt TM for Digital Micrograph from HREM Research, Japan (c.f. manuscript Fig. 1). For S Co(10), we can clearly see the evidence for Bloch DWs in Figure S7c.
In contrast, for Fe(0)/Co(10), due to the finite tilt angle, a part of the out-of-plane magnetization of the textures are projected in-plane and is visible in Figure S7d. Crucially, no evidence of Bloch DWs is detectable for Fe(0)/Co (10).
Domain compressibility investigated with LTEM. In manuscript Fig. 2, we present MTXM results as direct experimental evidence when investigating domain compressibility as a function of κ. The same trend could be observed indirectly with LTEM. To this end, LTEM imaging was performed along a hysteresis loop, as shown in Figure S8a-d. It is visually apparent that at low fields of −20 mT, the average domain width is much greater for Fe(0)/Co(10) ( Figure S8a) than that of Fe(2)/Co(8) ( Figure S8c). However, at higher fields, their domain widths are comparable ( Figure S8b,d). This visual observation could be placed on quantitative standing if we employ the previously established algorithmic linecut analysis, but this time as a function of the magnetic field, as shown in Figure S8e,f (c.f. MTXM in manuscript Fig. 2i). As the linecuts are centered on the intensity peaks, the domain width is simply the separation between the intensity peak (bright) and trough (dark) regions. For Fe(0)/Co(10) at low fields, the trough is separated ∼200 nm from the peak and is heavily smeared. This is because DWs are well separated (i.e. domains are wide) with weak spatial correlation. With increasing field, the trough moves towards the peak (i.e. domain width shrinks) and becomes more correlated. This is consistent with the highly compressible behavior seen in MTXM data. In contrast, for Fe(2)/Co(8), the peaks and troughs are always adjacent and have strong spatial correlation at all fields. This directly supports our claim that the domains are incompressible for

S4. Micromagnetic Simulations
Micromagnetic simulations were performed with MuMax3, using the magnetic parameters listed in Table S1 8 .
Layer-dependent chirality. In the manuscript, we discuss the lack of evidence for hybrid chirality in Fe(0)/Co(10), Fe(2)/Co(8) and Fe(3)/Co(7) based on our LTEM data. Our micromagnetic simulations -which will be elaborated here -are consistent with our experiments in suggesting that hybrid chirality, if present at all, is very limited in these three samples.
The zero-field micromagnetic simulations were analyzed and DWs regions for each layer were isolated for analysis. Figure S9a-d shows the OP magnetization (grayscale) and IP orientation of the DW magnetization relative However, upon considering all layers as a whole, both Bloch and Néel chiralities average to zero. Therefore, we label this sample as achiral 11,12 . In comparison, the other three samples -Fe(0)/Co(10), Fe(2)/Co (8) and Fe(3)/Co (7) -with D est > 0 are strongly chiral, i.e. have fixed handedness. Even if a Bloch center is present in Fe(0)/Co(10) (D est = 1.3 mJ/m 2 ), it might only exist in one or two layers which is difficult to observe experimentally. Moreover, we show that Fe(0)/Co(10) is at the threshold of hybrid chirality since the Bloch center vanishes if we consider a small but finite interlayer RKKY coupling -amounting to 20% of A est -as shown in Figure S9f  Simulation of hysteresis loops. Hysteresis loops are simulated by time-stepping in the presence of a sweeping applied magnetic field 8 . The rate of field sweep is approximately 10 6 T/s due to computational constraints, and so an entire hysteresis loop simulation is swept in about 600 ns. In order to ensure that simulated magnetization configurations may cross energy barriers at rates corresponding to conventional magnetometry experiments despite these constraints, the simulation temperature needs to be correspondingly elevated (e.g. to 850 -900 K). A typical hysteresis loop simulation result is shown in Figure S10 for direct comparison with experiments. Notwithstanding quantitative discrepancies of ∼20% in saturation field and magnetization, the simulated hysteresis loops fully reproduce the key experimental features, such as the sheared shape and various kinks.  Domain mapping. It is convenient to visualize skyrmion formation from magnetic stripes in simulations with a "field evolution map" of a typical magnetic domain. This is created by stacking the 2D footprint of a typical domain across varying fields into a 3D object, where the third axis is the magnetic field. This process is illustrated in Figure S11, where the chosen domain is in color. This 3D object is then rotated such that the field axis is horizontal, as shown in manuscript Fig. 4a,b. The mapping of the chosen domain across image slices in the magnetic field sweep is established using their spatial overlap -since the domains do not move appreciably within one field step.  Table S1 into model of isolated skyrmions 17 . This is compared with experimentally measured compressibility from manuscript Fig. 2k. The trends of both metrics are in close agreement -a sharp divide separates bubble skyrmions ( S Co(10), Fe(0)/Co(10)) from compact skyrmions (Fe(2)/Co(8), Fe(3)/Co (7)).
We also include two supplementary videos showing the evolution of the magnetic textures as a function of applied field that is used in the construction of manuscript Fig. 4(a,b).
SV1. This video shows the simulated evolution of magnetic textures of Fe(2)/Co(8) over a 2 µm field-of-view from 72 mT to 132 mT. Highlighted domains denote the representative evolution of stripes to skyrmions for Fe(2)/Co (8), and are horizontally stacked across fields to form manuscript Fig. 4(a).
SV2. This video shows the simulated evolution of magnetic textures of Fe(3)/Co(7) over a 2 µm field-of-view from 80 mT to 170 mT. Highlighted domains denote the representative evolution of stripes to skyrmions for Fe(3)/Co (7), and are horizontally stacked across fields to form manuscript Fig. 4(b).

S5. Geodesic Nudged Elastic Band (GNEB) simulations
GNEB atomistic calculations were performed with the Fidimag package 14 . A single layer of magnetic spins arranged in a two-dimensional square lattice with a cell size of 1 nm × 1 nm was used for all the GNEB simulations.
Energy barrier. GNEB is a well-established method to calculate the energy barrier for a transition between two fixed metastable magnetic configurations, under the constraint of fixed magnetic moments magnitude 15 . Since GNEB is an atomistic calculation, the energies of the configurations are calculated in accordance with the atomistic Heisenberg Hamiltonian, taking into account exchange, DMI, uniaxial anisotropy and magnetostatic interactions 16 .
To estimate the energy barriers of the stripe-skyrmion fission process, the initial metastable configuration was chosen to be single stripe domain relaxed at −0.04 mT, and the final configuration to be two metastable skyrmions

S6. Compressibility
Bubble skyrmions vs compact skyrmions. Here, we compare two metrics that can be used to differentiate bubble skyrmions and compact skyrmions. In Figure S12, we show the compressibility of stripes reproduced from manuscript Fig. 2k and minimum diameter of isolated skyrmions 17 for the four samples in this work. Both met-  (7) (compact skyrmions).
Correlation with K eff . One might notice a correlation between domain compressibility, ⟨dW/dH⟩ and K eff . We believe that this correlation arises from the fact that κ = πD/4 √ AK eff , and therefore K eff is a component of κ. Here, we show with micromagnetic simulations that K eff cannot fully account for the domain compressibility trends.
We study a fictitious sample, labeled as Fe(2)/Co(8)+, which has identical magnetic parameters as Fe(2)/Co(8)including K eff -except a direct exchange, A est , that is 20% larger, which results in a smaller κ. We simulate the hysteresis loop of Fe(2)/Co(8)+, extract the domain width and compute the ⟨dW/dH⟩ using the same protocol as detailed in the manuscript. The comparison between Fe(2)/Co (8) and Fe(2)/Co(8)+ is shown in Figure S13. As can be seen, ⟨dW/dH⟩ can vary considerably even if K eff is kept constant across samples.